function [fun, y0, A, F] = stelsel(a, b, c)
    %a = 1;
    %b = 1;
    %c = 1;
    A = zeros(3, 3);
    A(1,2) = 1;
    A(2,3) = 1;
    A(3,1) = -a * (b^2 + c^2);
    A(3,2) = -(b^2 + c^2 + 2 * a * b);
    A(3,3) = -(2*b + a);
    F = zeros(3,1);
    F(3) = a * (b^2 + c^2);
    
    fun = @(t, y) A * y + F;
    y0 = [2; c - a; a^2 - 2 * b * c];
end
